0. Setup
Exploratory Statistics
Define some plotting variables:
1. D1 Growth
mu_y = b3 + b4(diameter) + b5(height)
These models are actually doing a terrible job of predicting future diameter and height!
Which makes sense, because this model is supposed to be predicting growth of individuals that REMAIN in the seedling size domain and DO NOT recruit into the D2 size range.
2. Corrected to only include individuals within D1
–g1_1c Diameter only
–g1_5c Height only

0. Setup

Load the data

#Load demographic data
demog<-read.csv("demography_15_clean.csv", head=T)
library(AICcmodavg)
library(fields)
library(viridis)
library(scales)
my.palette <- viridis(9)

This is a csv file where each row is a an individual.

The column names are :

Obs : Observation number (each plant should have a different observation number; these were assigned post data entry when we discovered that some of the numbers for individuals were reused between sites
Site: Which site (Big Cypress, Cape Canaveral, Chekika, Forty Pierce, Punta Gorda, Wild Turkey)
Plant ID: The tag number of an individual (each individual in a site has a different number, but numbers were reused across sites)
Location: SeedlingPlot (individual located in a seedling plot) or NA (individual not in seedling plot)
Date_i: The date the census took place for censuses i = 1:7
Status_i: One of four categories for censuses i=1:7:
- alive: individual was alive when censused
- firstdead: the individual was recorded as dead at this census and was alive at the last census
- missing: the individual was not located at this census
- prestudy: the individual had not been tagged yet
- stilldead: the individual was already marked dead at a previous census
Death_Cause_i: One of three categories for censuses i = 1:7:
- NA: for individuals that died, there was no reason to suspect their deaths were due to management actions
- sprayed: individual death attributed to being sprayed with herbicide (these individuals should be excluded)
- fire: individual death attributed to fire (these individuals should be excluded from analysis)
Rep_stat_i: One of 4 categories for censuses i = 1:7:
- NA: for individuals whose reproductive status is unknown at census i
- female: for individuals who have produced fruit at any previous census up to an including census i
- male: for individuals who have not produced fruit at any previous census (may have observed male flowers)
- pre-reproductive: for individuals who have not produced fruit at any previous census
Diam_i: Diameter (in mm) of individual at census i
Height_i: Height (in cm) of individual at census i
Max_i: Maximum canopy width in cm of individual at census i
Min_i: Minimum canopy width in cm of individual at census i
Multistem: For each individual takes on one of two values:
- NA: individuals that are not multistem
- FLAG: Individuals that are multistem

Stack data

#### Restructure the data in a few different ways: 


Dates <- c(rep(2009, 1512), rep(2010, 1512), rep(2011, 1512), rep(2012, 1512), rep(2013, 1512), rep(2014, 1512))



Death_Cause<-c(demog$Death_Cause_2, demog$Death_Cause_3, demog$Death_Cause_4, demog$Death_Cause_5, demog$Death_Cause_6, demog$Death_Cause_7)

Status<-c(demog$Status_1, demog$Status_2, demog$Status_3, demog$Status_4, demog$Status_5, demog$Status_6, demog$Status_7)

Diams<-cbind(demog$Diam_1, demog$Diam_2, demog$Diam_3, demog$Diam_4, demog$Diam_5, demog$Diam_6, demog$Diam_7) 

row.names(Diams)<-demog$Plant_ID

Heights<-cbind(demog$Height_1, demog$Height_2, demog$Height_3, demog$Height_4, demog$Height_5, demog$Height_6, demog$Height_7)

row.names(Heights)<-demog$Plant_ID

Mins<-cbind(demog$Min_1, demog$Min_2, demog$Min_3, demog$Min_4, demog$Min_5, demog$Min_6, demog$Min_7)

row.names(Mins)<-demog$Plant_ID

Maxs<-cbind(demog$Max_1, demog$Max_2, demog$Max_3, demog$Max_4, demog$Max_5, demog$Max_6, demog$Max_7)

row.names(Maxs)<-demog$Plant_ID

Diams_pooled<-rbind( cbind(Diams[,1], Diams[,2]), cbind(Diams[,2], Diams[,3]), cbind(Diams[, 3], Diams[,4]), cbind(Diams[,4], Diams[,5]), cbind(Diams[,5], Diams[,6]), cbind(Diams[, 6], Diams[, 7]))

Diams_pooled<-data.frame(Diams_pooled[,1], Diams_pooled[,2])

names(Diams_pooled)<-c("t", "t_plus_one")

Sites_pooled<-rep(demog$Site, 6)

ID_pooled<-rep(demog$Plant_ID, 6)

Location<-rep(demog$Location, 6)

Genetic_type_pooled<-rep(0, length(Sites_pooled))

#code which sites are which biotype
for (i in 1:length(Genetic_type_pooled)) { 
    if (Sites_pooled[i]=="Big Cypress") Genetic_type_pooled[i]<-"hybrid"
    else if (Sites_pooled[i]=="Cape Canaveral") Genetic_type_pooled[i]<-"hybrid"
    else if (Sites_pooled[i]=="Chekika") Genetic_type_pooled[i]<-"eastern"
    else if (Sites_pooled[i]=="Fort Pierce") Genetic_type_pooled[i]<-"eastern"
    else if (Sites_pooled[i]=="Punta Gorda") Genetic_type_pooled[i]<-"western"
    else Genetic_type_pooled[i]<-"western"
    
}



Heights_pooled<-rbind( cbind(Heights[,1], Heights[,2]), cbind(Heights[,2], Heights[,3]), cbind(Heights[, 3], Heights[,4]), cbind(Heights[,4], Heights[,5]), cbind(Heights[,5], Heights[,6]), cbind(Heights[, 6], Heights[, 7]))

Heights_pooled<-data.frame(Heights_pooled[,1], Heights_pooled[,2])
names(Heights_pooled)<-c("t", "t_plus_one")



#Finish for maxs and mins 

Maxs_pooled<-rbind(cbind(Maxs[,1], Maxs[,2]), cbind(Maxs[,2], Maxs[,3]), 
    cbind(Maxs[,3], Maxs[,4]), cbind(Maxs[,4], Maxs[,5]), 
    cbind(Maxs[,5], Maxs[,6]), cbind(Maxs[,6], Maxs[,7]))
    
Maxs_pooled<-data.frame(Maxs_pooled[,1], Maxs_pooled[,2])
names(Maxs_pooled)<-c("t", "t_plus_one")



Mins_pooled<-rbind( cbind(Mins[,1], Mins[,2]), cbind(Mins[,2], Mins[,3]), 
    cbind(Mins[,3], Mins[,4]), cbind(Mins[,4], Mins[,5]), 
    cbind(Mins[, 5], Mins[,6]), cbind(Mins[,6], Mins[,7]))
    
Mins_pooled<-data.frame(Mins_pooled[,1], Mins_pooled[,2])
names(Mins_pooled)<-c("t", "t_plus_one")




###Use status matrix to get survival matrix
statmat<-as.matrix(demog[, c(6, 14, 22, 30, 38, 46, 54)])
survmat<-statmat
survmat<- ifelse(survmat=="prestudy"|survmat=="missing"|survmat=="still_dead", NA, survmat)
survmat<- ifelse(survmat =="tagged"|survmat=="alive", 1, 0)

survmat_pooled<-rbind( cbind(survmat[,1], survmat[,2]), cbind(survmat[,2], survmat[,3]), cbind(survmat[, 3], survmat[,4]), cbind(survmat[,4], survmat[,5]), cbind(survmat[,5], survmat[,6]), cbind(survmat[, 6], survmat[, 7]))
names(survmat_pooled)<-c("surv_t", "surv_t_plus_one")


survmat_pooled_bysize<-cbind(Diams_pooled[, 1], survmat_pooled[,2])
names(survmat_pooled_bysize)<-c("size_t", "surv_t_plus_one")



#Use reproduction status matrix to get reproduction matrix
repstatmat<-as.matrix(demog[, c(8, 16, 24, 32, 40, 48, 56)])
repmat<-repstatmat
repmat<- ifelse(repmat =="female", 1, 0)
repmat_pooled<-rbind( cbind(repmat[,1], repmat[,2]), cbind(repmat[,2], repmat[,3]), cbind(repmat[, 3], repmat[,4]), cbind(repmat[,4], repmat[,5]), cbind(repmat[,5], repmat[,6]), cbind(repmat[, 6], repmat[, 7]))


pooled<-cbind(ID_pooled, Sites_pooled, Location, Death_Cause, Genetic_type_pooled, Mins_pooled, Maxs_pooled, Diams_pooled, Heights_pooled, repmat_pooled, survmat_pooled, Dates)
names(pooled)<-c("ID", "Site","Location", "Death_Cause", "Genetic_type", "Min_t", "Min_tplus1", "Max_t", "Max_tplus1", "Diameter_t", "Diameter_tplus1", "Height_t", "Height_tplus1", "Rep_t", "Rep_tplus1", "Surv_t", "Surv_tplus1", "Year")

Investigate relationships among possible state variables

predictors<-cbind(pooled$Min_t, pooled$Diameter_t, pooled$Height_t, pooled$Max_t)

# What is the correlation amongst the predictor variables? 
correlation<-cor(predictors[, 1:4], use="pairwise.complete.obs")
colnames(correlation)<-c("Min", "Diameter", "Height", "Max")
rownames(correlation)<-c("Min", "Diameter", "Height", "Max")

correlation<-round(correlation, 2)
print(correlation)

##           Min Diameter Height  Max
## Min      1.00     0.70   0.85 0.94
## Diameter 0.70     1.00   0.58 0.68
## Height   0.85     0.58   1.00 0.87
## Max      0.94     0.68   0.87 1.00

#Over the entire dataset, there is the lowest correlation between diameter and height

Exclude outliers

#Remove individual deaths that were sprayed or fire or grinding
restrict<-subset(pooled, is.na(Death_Cause))

#Remove outliers (very large diameters)
restrict2<-subset(restrict, Diameter_t<800)

#Remove an individual at Big Cypress that went from diameter 1.4 to diameter 108 back to diameter 2.1 
restrict2 <- subset(restrict2, restrict2$ID != 5)

#Questioning why we did this? does this lead to model-eviction? 
#Range of diameter: [0, 800]
0.9*max(restrict2$Diameter_t) #700

## [1] 625.5

#Range of height:  [0, 800]
0.9*max(restrict2$Height_t, na.rm=T) #800

## [1] 729

Divide dataset into size domains (seedlings and larges)

seedlings<-subset(restrict2, restrict2$Diameter_t<1.6)
seedlings<-subset(seedlings, seedlings$Height_t<16)
larges<-subset(restrict2, restrict2$Diameter_t>1.6)
larges<-subset(larges, larges$Height_t>16)

x1seq_larges<-seq(1.6, 800, length.out=50)
x2seq_larges<-seq(16, 800, length.out=50)

x1seq_seedlings<-seq(0, 1.6, length.out=50)
x2seq_seedlings<-seq(0, 16, length.out=50)

grad_status<-rep(NA, length(seedlings$ID))
for (i in 1:length(seedlings$ID)) {
  if (is.na(seedlings$Diameter_tplus1[i])) {
    grad_status[i]<-NA
  }
  else if (is.na(seedlings$Height_tplus1[i])) {
    grad_status[i]<-NA
  }
  else if (seedlings$Diameter_tplus1[i]>1.6& seedlings$Height_tplus1[i]>16) { grad_status[i]<-1}
  else {grad_status[i]<-0}
}

seedlings<-cbind(seedlings, grad_status)

seedlings_g<-subset(seedlings, seedlings$Diameter_tplus1<1.6)
seedlings_g<-subset(seedlings_g, seedlings_g$Height_tplus1<16)

Exploratory Statistics

There are 1512 individuals

Big Cypress	Cape Canaveral	Chekika	Fort Pierce	Punta Gorda	Wild Turkey
308	226	460	166	175	177

Eastern	Hybrid	Western
626	534	352

summary(seedlings$Site)

##    Big Cypress Cape Canaveral        Chekika    Fort Pierce    Punta Gorda 
##            209            165            389             76             58 
##    Wild Turkey 
##             60

summary(seedlings$Genetic_type)

## eastern  hybrid western 
##     465     374     118

summary(as.factor(seedlings$Year))

## 2009 2010 2011 2012 2013 2014 
##  119  194   59  351  165   69

summary(larges$Site)

##    Big Cypress Cape Canaveral        Chekika    Fort Pierce    Punta Gorda 
##            680            483            546            542            422 
##    Wild Turkey 
##            477

summary(larges$Genetic_type)

## eastern  hybrid western 
##    1088    1163     899

summary(as.factor(larges$Year))

## 2009 2010 2011 2012 2013 2014 
##  299  631  644  571  518  487

The number of deaths (0) by biotype:

table(seedlings$Surv_tplus1, seedlings$Genetic_type)

##    
##     eastern hybrid western
##   0     304    165      67
##   1     146    168      43

table(larges$Surv_tplus1, larges$Genetic_type)

##    
##     eastern hybrid western
##   0      82     33      66
##   1     983    921     822

The number of deaths(0) by site:

table(seedlings$Surv_tplus1, seedlings$Site)

##    
##     Big Cypress Cape Canaveral Chekika Fort Pierce Punta Gorda Wild Turkey
##   0         119             46     277          27          31          36
##   1          71             97     104          42          20          23

table(larges$Surv_tplus1, larges$Site)

##    
##     Big Cypress Cape Canaveral Chekika Fort Pierce Punta Gorda Wild Turkey
##   0          19             14      51          31          20          46
##   1         549            372     478         505         397         425

table(seedlings$Surv_tplus1, seedlings$Year)

##    
##     2009 2010 2011 2012 2013 2014
##   0   24   81   36  264  120   11
##   1   95  104   23   86   45    4

table(larges$Surv_tplus1, larges$Year)

##    
##     2009 2010 2011 2012 2013 2014
##   0    0   31   51   49   28   22
##   1  299  595  581  517  489  245

Define some plotting variables:

x1seq_seedlings<-seq(0, 1.6, length.out=50)

x2seq_seedlings<-seq(0, 16, length.out=50)
x1seq_larges <- seq(0, 800, length.out = 50 )
x2seq_larges <- seq (0, 800, length.out = 50)

invlogit<-function(x){exp(x)/(1+exp(x))}

# Function to do an image plot of a matrix in the usual orientation, A(1,1) at top left  
matrix.image=function(x,y,A,col=topo.colors(200),...) {
    nx=length(x); ny=length(y); 
    x1=c(1.5*x[1]-0.5*x[2],1.5*x[nx]-0.5*x[nx-1]); 
    y1=c(1.5*y[1]-0.5*y[2],1.5*y[ny]-0.5*y[ny-1]); 
    image.plot(list(x=x,y=y,z=t(A)),xlim=x1,ylim=rev(y1),col=col,bty="u",...);  
}

1. D1 Growth

Choosing which state variables to use and whether to include terms for biotype, site or year (for the year effect, only including data from years where all sites were sampled, so AIC cannot be compared with other models)

Diameter-only models:
- m1: logit(m) = b0 + b1(diameter)
- m2: logit(m) = b0 + b1(diameter) + epsilon(biotype{E,H,W})
- m3: logit(m) = b0 + b1(diameter) + epsilon(site{BC, CC, C, FP, PG, WT})
- m4: logit(m) = b0 + b1(diameter) + epsilon(year{2010, 2011, 2012, 2013, 2014}
Height-only models
- m5: logit(m) = b0 + b1(height)
- m6: logit(m) = b0 + b1(height) + epsilon(biotype{E,H,W})
- m7: logit(m) = b0 + b1(height) + epsilon(site{BC, CC, C, FP, PG, WT})
- m8: logit(m) = b0 + b1(height) + epsilon(year{2010, 2011, 2012, 2013, 2014}
Diameter + Height models
- m9: logit(m) = b0 + b1(diameter) + b2(height)
- m10: logit(m) = b0 + b1(diameter) + b2(height) + epsilon(biotype{E,H,W})
- m11: logit(m) = b0 + b1(diameter) + b2(height) + epsilon(site{BC, CC, C, FP, PG,WT})
- m12: loigt(m) = b0 + b1(diameter) + b2(height) + epsilon(year{2010, 2011, 2012, 2013, 2014})

-Diameter-only models

– g1_1: mu_x= b0 + b1(diameter) —–

g1_1 <- lm(seedlings_g$Diameter_tplus1 ~ seedlings_g$Diameter_t)
summary(g1_1)

## 
## Call:
## lm(formula = seedlings_g$Diameter_tplus1 ~ seedlings_g$Diameter_t)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83841 -0.12663 -0.01484  0.10569  0.54105 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.56769    0.04895   11.60   <2e-16 ***
## seedlings_g$Diameter_t  0.55894    0.05317   10.51   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1927 on 243 degrees of freedom
## Multiple R-squared:  0.3126, Adjusted R-squared:  0.3098 
## F-statistic: 110.5 on 1 and 243 DF,  p-value: < 2.2e-16

#Variance in Growth
sd_x <- summary(g1_1)$sigma


preds2 <- g1_1$coefficients[1] + g1_1$coefficients[2]*x1seq_seedlings + 2*sd_x
preds3 <- g1_1$coefficients[1] + g1_1$coefficients[2]*x1seq_seedlings - 2*sd_x
preds4 <- g1_1$coefficients[1] + g1_1$coefficients[2]*x1seq_seedlings + 1*sd_x
preds5 <- g1_1$coefficients[1] + g1_1$coefficients[2]*x1seq_seedlings - 1*sd_x



plot(x1seq_seedlings, g1_1$coeff[1] + g1_1$coeff[2]*x1seq_seedlings, type='l' , xlab = "Diameter at t (mm)", ylab="Diameter at t+1 (mm)")
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings_g$Diameter_t, seedlings_g$Diameter_tplus1)
abline(g1_1)

–g1_2: mu_x = b0 + b1(diameter) + epsilon(biotype{E,H,W})

Hybrid is not significantly different from Eastern

g1_2 <- lm(seedlings_g$Diameter_tplus1 ~ seedlings_g$Diameter_t + seedlings_g$Genetic_type)
summary(g1_2)

## 
## Call:
## lm(formula = seedlings_g$Diameter_tplus1 ~ seedlings_g$Diameter_t + 
##     seedlings_g$Genetic_type)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.82788 -0.11841 -0.00767  0.10472  0.46395 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     0.563404   0.048832  11.537   <2e-16 ***
## seedlings_g$Diameter_t          0.553729   0.054492  10.162   <2e-16 ***
## seedlings_g$Genetic_typehybrid  0.001281   0.026564   0.048   0.9616    
## seedlings_g$Genetic_typewestern 0.085037   0.043626   1.949   0.0524 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1919 on 241 degrees of freedom
## Multiple R-squared:  0.3243, Adjusted R-squared:  0.3159 
## F-statistic: 38.56 on 3 and 241 DF,  p-value: < 2.2e-16

#Variance in Growth
sd_x <- summary(g1_2)$sigma

par(mfrow=c(1,2))
#Eastern/Hybrid
preds2 <- g1_2$coefficients[1] + g1_2$coefficients[2]*x1seq_seedlings  + 2*sd_x
preds3 <- g1_2$coefficients[1] + g1_2$coefficients[2]*x1seq_seedlings - 2*sd_x
preds4 <- g1_2$coefficients[1] + g1_2$coefficients[2]*x1seq_seedlings + 1*sd_x
preds5 <- g1_2$coefficients[1] + g1_2$coefficients[2]*x1seq_seedlings - 1*sd_x

plot(x1seq_seedlings, g1_2$coeff[1] + g1_2$coeff[2]*x1seq_seedlings, type='l' , xlab = "Diameter at t (mm)", ylab="Diameter at t+1 (mm)")
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings_g$Diameter_t[seedlings_g$Genetic_type == "eastern" | seedlings_g$Genetic_type == "hybrid"], seedlings_g$Diameter_tplus1[seedlings_g$Genetic_type == "eastern" | seedlings_g$Genetic_type == "hybrid"])
lines(x1seq_seedlings, g1_2$coeff[1] + g1_2$coeff[2]*x1seq_seedlings)

#Western
preds2 <- g1_2$coefficients[1] + g1_2$coefficients[2]*x1seq_seedlings + g1_2$coefficients[4]   + 2*sd_x
preds3 <- g1_2$coefficients[1] + g1_2$coefficients[2]*x1seq_seedlings +   
  g1_2$coefficients[4]- 2*sd_x
preds4 <- g1_2$coefficients[1] + g1_2$coefficients[2]*x1seq_seedlings + 
  g1_2$coefficients[4]+ 1*sd_x
preds5 <- g1_2$coefficients[1] + g1_2$coefficients[2]*x1seq_seedlings + 
  g1_2$coefficients[4]- 1*sd_x

plot(x1seq_seedlings, g1_2$coeff[1] + g1_2$coeff[2]*x1seq_seedlings + g1_2$coefficients[4], 
     type='l' , xlab = "Diameter at t (mm)", ylab="Diameter at t+1 (mm)")
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds3), preds2), col = 'grey90', 
        border = NA)
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds5), preds4), col = 'grey80', 
        border = NA)
points(seedlings_g$Diameter_t[seedlings_g$Genetic_type == "western"], seedlings_g$Diameter_tplus1[seedlings_g$Genetic_type == "western"])
lines(x1seq_seedlings, g1_2$coeff[1] + g1_2$coeff[2]*x1seq_seedlings + g1_2$coefficients[4])

–g1_3: mu_x= b0 + b1(diameter) + epsilon(biotype{BC, CC, C, FP, PG, WT})

Big Cypress, Cape Canaveral and Chekika are not different

g1_3 <- lm(seedlings_g$Diameter_tplus1 ~ seedlings_g$Diameter_t + seedlings_g$Site)
summary(g1_3)

## 
## Call:
## lm(formula = seedlings_g$Diameter_tplus1 ~ seedlings_g$Diameter_t + 
##     seedlings_g$Site)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.80112 -0.12797 -0.01872  0.08919  0.50085 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     0.52201    0.05603   9.316  < 2e-16 ***
## seedlings_g$Diameter_t          0.54628    0.05550   9.844  < 2e-16 ***
## seedlings_g$SiteCape Canaveral  0.07868    0.03741   2.103  0.03651 *  
## seedlings_g$SiteChekika         0.02356    0.03637   0.648  0.51769    
## seedlings_g$SiteFort Pierce     0.10640    0.04526   2.351  0.01954 *  
## seedlings_g$SitePunta Gorda     0.17884    0.06432   2.781  0.00586 ** 
## seedlings_g$SiteWild Turkey     0.09474    0.06005   1.578  0.11598    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1892 on 238 degrees of freedom
## Multiple R-squared:  0.3515, Adjusted R-squared:  0.3351 
## F-statistic:  21.5 on 6 and 238 DF,  p-value: < 2.2e-16

#Variance in Growth
sd_x <- summary(g1_3)$sigma

par(mfrow=c(2,2))
#Big Cypress, Cape Canaveral and Chekika 
preds2 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings  + 2*sd_x
preds3 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings - 2*sd_x
preds4 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 1*sd_x
preds5 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings - 1*sd_x

plot(x1seq_seedlings, g1_3$coeff[1] + g1_3$coeff[2]*x1seq_seedlings, type='l' , xlab = "Diameter at t (mm)", ylab="Diameter at t+1 (mm)", main = "BC, CC, C")
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings_g$Diameter_t[seedlings_g$Site == "Big Cypress" | seedlings_g$Site == "Cape Canaveral" | seedlings_g$Site == "Chekika"], seedlings_g$Diameter_tplus1[seedlings_g$Site == "Big Cypress" | seedlings_g$Site == "Cape Canaveral" | seedlings_g$Site == "Chekika"])
lines(x1seq_seedlings, g1_3$coeff[1] + g1_3$coeff[2]*x1seq_seedlings)

#Fort Pierce
preds2 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 
  g1_3$coefficients[5]  + 2*sd_x
preds3 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings +   
  g1_3$coefficients[5]- 2*sd_x
preds4 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 
  g1_3$coefficients[5]+ 1*sd_x
preds5 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 
  g1_3$coefficients[5]- 1*sd_x

plot(x1seq_seedlings, g1_3$coeff[1] + g1_3$coeff[2]*x1seq_seedlings + g1_3$coefficients[5], 
     type='l' , xlab = "Diameter at t (mm)", ylab="Diameter at t+1 (mm)", main = "FP")
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds3), preds2), col = 'grey90', 
        border = NA)
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds5), preds4), col = 'grey80', 
        border = NA)
points(seedlings_g$Diameter_t[seedlings_g$Site == "Fort Pierce"], seedlings_g$Diameter_tplus1[seedlings_g$Site == "Fort Pierce"])
lines(x1seq_seedlings, g1_3$coeff[1] + g1_3$coeff[2]*x1seq_seedlings + g1_3$coefficients[5])

#Punta Gorda
preds2 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 
  g1_3$coefficients[6]   + 2*sd_x
preds3 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings +   
  g1_3$coefficients[6]- 2*sd_x
preds4 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 
  g1_3$coefficients[6]+ 1*sd_x
preds5 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 
  g1_3$coefficients[6]- 1*sd_x

plot(x1seq_seedlings, g1_3$coeff[1] + g1_3$coeff[2]*x1seq_seedlings + g1_3$coefficients[6], 
     type='l' , xlab = "Diameter at t (mm)", ylab="Diameter at t+1 (mm)", main = "PG")
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds3), preds2), col = 'grey90', 
        border = NA)
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds5), preds4), col = 'grey80', 
        border = NA)
points(seedlings_g$Diameter_t[seedlings_g$Site == "Punta Gorda"], seedlings_g$Diameter_tplus1[seedlings_g$Site == "Punta Gorda"])
lines(x1seq_seedlings, g1_3$coeff[1] + g1_3$coeff[2]*x1seq_seedlings + g1_3$coefficients[6])

#Wild Turkey
preds2 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 
  g1_3$coefficients[7]   + 2*sd_x
preds3 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings +   
  g1_3$coefficients[7]- 2*sd_x
preds4 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 
  g1_3$coefficients[7]+ 1*sd_x
preds5 <- g1_3$coefficients[1] + g1_3$coefficients[2]*x1seq_seedlings + 
  g1_3$coefficients[7]- 1*sd_x

plot(x1seq_seedlings, g1_3$coeff[1] + g1_3$coeff[2]*x1seq_seedlings + g1_3$coefficients[7], 
     type='l' , xlab = "Diameter at t (mm)", ylab="Diameter at t+1 (mm)", main = "WT")
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds3), preds2), col = 'grey90', 
        border = NA)
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds5), preds4), col = 'grey80', 
        border = NA)
points(seedlings_g$Diameter_t[seedlings_g$Site == "Wild Turkey"], seedlings_g$Diameter_tplus1[seedlings_g$Site == "Wild Turkey"])
lines(x1seq_seedlings, g1_3$coeff[1] + g1_3$coeff[2]*x1seq_seedlings + g1_3$coefficients[7])

–g1_4: mu_x = b0 + b1(diameter) + epsilon(year{2010, 2011, 2012, 2013}

There is no difference in seedling growth among years.

seedlings2 <- subset(seedlings, seedlings$Year > 2009 & seedlings$Year <2014)
seedlings2$Year <- as.factor(seedlings2$Year)
g1_4 <- lm(seedlings2$Diameter_tplus1 ~ seedlings2$Diameter_t + seedlings2$Year)
summary(g1_4)

## 
## Call:
## lm(formula = seedlings2$Diameter_tplus1 ~ seedlings2$Diameter_t + 
##     seedlings2$Year)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.03692 -0.20882 -0.06532  0.11829  2.36308 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.64930    0.12599   5.154 5.23e-07 ***
## seedlings2$Diameter_t  0.90635    0.11696   7.749 2.40e-13 ***
## seedlings2$Year2011   -0.38434    0.09820  -3.914 0.000118 ***
## seedlings2$Year2012   -0.39033    0.06235  -6.261 1.69e-09 ***
## seedlings2$Year2013   -0.35919    0.07419  -4.842 2.27e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4157 on 247 degrees of freedom
##   (517 observations deleted due to missingness)
## Multiple R-squared:  0.3241, Adjusted R-squared:  0.3132 
## F-statistic: 29.61 on 4 and 247 DF,  p-value: < 2.2e-16

#Variance in Growth
sd_x <- summary(g1_4)$sigma

- Height-only models

–g1_5: mu_y = b0 + b1(height)

  g1_5 <- lm(seedlings$Height_tplus1 ~ seedlings$Height_t)
summary(g1_5)

## 
## Call:
## lm(formula = seedlings$Height_tplus1 ~ seedlings$Height_t)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.6293  -2.8908  -0.8908   1.9131  28.3053 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          3.5444     0.9325   3.801  0.00017 ***
## seedlings$Height_t   0.9346     0.0947   9.869  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.169 on 352 degrees of freedom
##   (603 observations deleted due to missingness)
## Multiple R-squared:  0.2167, Adjusted R-squared:  0.2145 
## F-statistic:  97.4 on 1 and 352 DF,  p-value: < 2.2e-16

#Variance in Growth
sd_x <- summary(g1_5)$sigma


preds2 <- g1_5$coefficients[1] + g1_5$coefficients[2]*x2seq_seedlings + 2*sd_x
preds3 <- g1_5$coefficients[1] + g1_5$coefficients[2]*x2seq_seedlings - 2*sd_x
preds4 <- g1_5$coefficients[1] + g1_5$coefficients[2]*x2seq_seedlings + 1*sd_x
preds5 <- g1_5$coefficients[1] + g1_5$coefficients[2]*x2seq_seedlings - 1*sd_x



plot(x2seq_seedlings, g1_5$coeff[1] + g1_5$coeff[2]*x2seq_seedlings, type='l' , xlab = "Height at t (cm)", ylab="Height at t+1 (cm)")
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings$Height_t, seedlings$Height_tplus1)
abline(g1_5)

–g1_6: mu_y = b0 + b1(height) + epsilon(biotype{E,H,W})

Hybrid and Eastern are not different

g1_6 <- lm(seedlings$Height_tplus1 ~ seedlings$Height_t + seedlings$Genetic_type)
summary(g1_6)

## 
## Call:
## lm(formula = seedlings$Height_tplus1 ~ seedlings$Height_t + seedlings$Genetic_type)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.7532  -3.0298  -0.4268   1.8057  27.1747 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    2.85965    0.98190   2.912  0.00382 ** 
## seedlings$Height_t             0.92787    0.09425   9.845  < 2e-16 ***
## seedlings$Genetic_typehybrid   1.05599    0.58834   1.795  0.07354 .  
## seedlings$Genetic_typewestern  1.90338    0.87428   2.177  0.03014 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.141 on 350 degrees of freedom
##   (603 observations deleted due to missingness)
## Multiple R-squared:  0.2297, Adjusted R-squared:  0.2231 
## F-statistic: 34.79 on 3 and 350 DF,  p-value: < 2.2e-16

#Variance in Growth
sd_x <- summary(g1_6)$sigma

par(mfrow=c(1,2))
#Eastern/Hybrid
preds2 <- g1_6$coefficients[1] + g1_6$coefficients[2]*x2seq_seedlings  + 2*sd_x
preds3 <- g1_6$coefficients[1] + g1_6$coefficients[2]*x2seq_seedlings - 2*sd_x
preds4 <- g1_6$coefficients[1] + g1_6$coefficients[2]*x2seq_seedlings + 1*sd_x
preds5 <- g1_6$coefficients[1] + g1_6$coefficients[2]*x2seq_seedlings - 1*sd_x

plot(x2seq_seedlings, g1_6$coeff[1] + g1_6$coeff[2]*x2seq_seedlings, type='l' , xlab = "Height at t (cm)", ylab="Height at t+1 (cm)")
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings$Height_t[seedlings$Genetic_type == "eastern" | seedlings$Genetic_type == "hybrid"], seedlings$Height_tplus1[seedlings$Genetic_type == "eastern" | seedlings$Genetic_type == "hybrid"])
lines(x2seq_seedlings, g1_6$coeff[1] + g1_6$coeff[2]*x2seq_seedlings)

#Western
preds2 <- g1_6$coefficients[1] + g1_6$coefficients[2]*x2seq_seedlings + g1_6$coefficients[4]   + 2*sd_x
preds3 <- g1_6$coefficients[1] + g1_6$coefficients[2]*x2seq_seedlings +   
  g1_6$coefficients[4]- 2*sd_x
preds4 <- g1_6$coefficients[1] + g1_6$coefficients[2]*x2seq_seedlings + 
  g1_6$coefficients[4]+ 1*sd_x
preds5 <- g1_6$coefficients[1] + g1_6$coefficients[2]*x2seq_seedlings + 
  g1_6$coefficients[4]- 1*sd_x

plot(x2seq_seedlings, g1_6$coeff[1] + g1_6$coeff[2]*x2seq_seedlings + g1_6$coefficients[4], 
     type='l' , xlab = "Height at t (cm)", ylab="Height at t+1 (cm)")
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds3), preds2), col = 'grey90', 
        border = NA)
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds5), preds4), col = 'grey80', 
        border = NA)
points(seedlings$Height_t[seedlings$Genetic_type == "western"], seedlings$Height_tplus1[seedlings$Genetic_type == "western"])
lines(x2seq_seedlings, g1_6$coeff[1] + g1_6$coeff[2]*x2seq_seedlings + g1_6$coefficients[4])

–g1_7: mu_y = b0 + b1(height) + epsilon(site{BC, CC, C, FP, PG, WT})

Wild Turkey is the only site distict from all the rest (but how much of it is just that Wild Turkey has a much smaller sample size than the other sites?)

 g1_7 <- lm(seedlings_g$Height_tplus1 ~ seedlings_g$Height_t + seedlings_g$Site)
summary(g1_7)

## 
## Call:
## lm(formula = seedlings_g$Height_tplus1 ~ seedlings_g$Height_t + 
##     seedlings_g$Site)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.6087 -1.8749  0.0247  1.6861  6.0366 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     6.93205    0.80003   8.665 7.05e-16 ***
## seedlings_g$Height_t            0.35974    0.06987   5.149 5.49e-07 ***
## seedlings_g$SiteCape Canaveral -0.35284    0.54283  -0.650   0.5163    
## seedlings_g$SiteChekika         0.06495    0.52274   0.124   0.9012    
## seedlings_g$SiteFort Pierce    -0.36466    0.64254  -0.568   0.5709    
## seedlings_g$SitePunta Gorda     0.14291    0.91552   0.156   0.8761    
## seedlings_g$SiteWild Turkey    -2.20629    0.85562  -2.579   0.0105 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.693 on 238 degrees of freedom
## Multiple R-squared:  0.1306, Adjusted R-squared:  0.1087 
## F-statistic: 5.958 on 6 and 238 DF,  p-value: 8.103e-06

#Variance in Growth
sd_x <- summary(g1_7)$sigma

par(mfrow=c(1,2))
#All sites except Wild Turkey 
preds2 <- g1_7$coefficients[1] + g1_7$coefficients[2]*x2seq_seedlings  + 2*sd_x
preds3 <- g1_7$coefficients[1] + g1_7$coefficients[2]*x2seq_seedlings - 2*sd_x
preds4 <- g1_7$coefficients[1] + g1_7$coefficients[2]*x2seq_seedlings + 1*sd_x
preds5 <- g1_7$coefficients[1] + g1_7$coefficients[2]*x2seq_seedlings - 1*sd_x

plot(x2seq_seedlings, g1_7$coeff[1] + g1_7$coeff[2]*x2seq_seedlings, type='l' , xlab = "Height at t (mm)", ylab="Height at t+1 (mm)", main = "BC, CC, C, FP, PG")
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings_g$Height_t[seedlings_g$Site != "Wild Turkey"],
       seedlings_g$Height_tplus1[seedlings_g$Site != "Wild Turkey"])
lines(x2seq_seedlings, g1_7$coeff[1] + g1_7$coeff[2]*x2seq_seedlings)

plot(x2seq_seedlings, g1_7$coeff[1] + g1_7$coeff[2]*x2seq_seedlings +g1_7$coeff[7], type='l' , xlab = "Height at t (mm)", ylab="Height at t+1 (mm)", main = "WT")
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings_g$Height_t[seedlings_g$Site == "Wild Turkey"],
       seedlings_g$Height_tplus1[seedlings_g$Site == "Wild Turkey"])
lines(x2seq_seedlings, g1_7$coeff[1] + g1_7$coeff[2]*x2seq_seedlings + g1_7$coeff[7])

–g1_8: mu_y = b0 + b1(height) + epsilon(year{2010, 2011, 2012, 2013}

2013 is different from the rest of the years, although once again it could just be the smaller sample size.

g1_8 <- lm(seedlings2$Height_tplus1 ~ seedlings2$Height_t + seedlings2$Year)
summary(g1_8)

## 
## Call:
## lm(formula = seedlings2$Height_tplus1 ~ seedlings2$Height_t + 
##     seedlings2$Year)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.1909  -3.3749  -0.6768   2.4152  26.7588 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           4.8952     1.3689   3.576 0.000419 ***
## seedlings2$Height_t   0.9497     0.1234   7.697 3.22e-13 ***
## seedlings2$Year2011  -2.1291     1.2963  -1.642 0.101752    
## seedlings2$Year2012  -1.0172     0.8852  -1.149 0.251607    
## seedlings2$Year2013  -2.0266     1.0033  -2.020 0.044463 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.64 on 250 degrees of freedom
##   (514 observations deleted due to missingness)
## Multiple R-squared:  0.2342, Adjusted R-squared:  0.222 
## F-statistic: 19.12 on 4 and 250 DF,  p-value: 9.86e-14

#Variance in Growth
sd_x <- summary(g1_8)$sigma

preds2 <- g1_8$coefficients[1] + g1_8$coefficients[2]*x2seq_seedlings  + 2*sd_x
preds3 <- g1_8$coefficients[1] + g1_8$coefficients[2]*x2seq_seedlings - 2*sd_x
preds4 <- g1_8$coefficients[1] + g1_8$coefficients[2]*x2seq_seedlings + 1*sd_x
preds5 <- g1_8$coefficients[1] + g1_8$coefficients[2]*x2seq_seedlings - 1*sd_x

par(mfrow=c(1,2))
plot(x2seq_seedlings, g1_8$coeff[1] + g1_8$coeff[2]*x2seq_seedlings, type='l' , xlab = "Height at t (mm)", ylab="Height at t+1 (mm)", main = "2010, 2011, 2012")
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings2$Height_t[seedlings2$Year != "2013"],
       seedlings2$Height_tplus1[seedlings2$Year != "2013"])
lines(x2seq_seedlings, g1_8$coeff[1] + g1_8$coeff[2]*x2seq_seedlings)

preds2 <- g1_8$coefficients[1] + g1_8$coefficients[2]*x2seq_seedlings + g1_8$coeff[5]  + 2*sd_x
preds3 <- g1_8$coefficients[1] + g1_8$coefficients[2]*x2seq_seedlings+ g1_8$coeff[5] - 2*sd_x
preds4 <- g1_8$coefficients[1] + g1_8$coefficients[2]*x2seq_seedlings + g1_8$coeff[5]+ 1*sd_x
preds5 <- g1_8$coefficients[1] + g1_8$coefficients[2]*x2seq_seedlings + g1_8$coeff[5]- 1*sd_x
plot(x2seq_seedlings, g1_8$coeff[1] + g1_8$coeff[2]*x2seq_seedlings + g1_8$coeff[5], type='l' , xlab = "Height at t (mm)", ylab="Height at t+1 (mm)", main = "2013")
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings2$Height_t[seedlings2$Year == "2013"],
       seedlings2$Height_tplus1[seedlings2$Year == "2013"])
lines(x2seq_seedlings, g1_8$coeff[1] + g1_8$coeff[2]*x2seq_seedlings + g1_8$coeff[5])

### - Diameter + Height models

–g1_9: mu_x = b0 + b1(diameter) + b2(height)

mu_y = b3 + b4(diameter) + b5(height)

These models are actually doing a terrible job of predicting future diameter and height!

Which makes sense, because this model is supposed to be predicting growth of individuals that REMAIN in the seedling size domain and DO NOT recruit into the D2 size range.

  g1_9 <-manova(cbind(seedlings$Diameter_tplus1, seedlings$Height_tplus1) ~ seedlings$Diameter_t+seedlings$Height_t)
summary(g1_9)

##                       Df  Pillai approx F num Df den Df    Pr(>F)    
## seedlings$Diameter_t   1 0.27583   66.085      2    347 < 2.2e-16 ***
## seedlings$Height_t     1 0.16005   33.059      2    347 7.215e-14 ***
## Residuals            348                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sd_x<-summary(lm(seedlings$Diameter_tplus1 ~ seedlings$Diameter_t + seedlings$Height_t))$sigma

#Grab parameters associated with future diameter
b0<-g1_9$coeff[1,1]
b1<-g1_9$coeff[2,1]
b2<-g1_9$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_9$coeff[1,2]
b4<-g1_9$coeff[2,2]
b5<-g1_9$coeff[3,2]
sd_y <-summary(lm(seedlings$Height_tplus1 ~ seedlings$Diameter_t + seedlings$Height_t))$sigma
 
 

z1<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b)))
z2<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

par(mfrow=c(1,2))
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings$Diameter_t, seedlings$Height_t, seedlings$Diameter_tplus1, pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)
#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings$Diameter_t, seedlings$Height_t, seedlings$Height_tplus1, pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

par(mfrow=c(1,1))
par(mar=c(5.1,4.1,4.1,6))
coeff_bigger=1
image.plot(x1seq_seedlings, x2seq_seedlings, z1, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)")
points(seedlings$Diameter_t, seedlings$Height_t, cex=seedlings$Diameter_tplus1*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )


my.palette <- viridis(9)
par(mar=c(5.1,4.1,4.1,6))
coeff_bigger=0.1
image.plot(x1seq_seedlings, x2seq_seedlings, z2, xlab = "Diameter (mm) ", ylab="Height (cm)", col=alpha(my.palette, 0.2), legend.lab = "Height at t+1 (cm)")
points(seedlings$Diameter_t, seedlings$Height_t, cex=seedlings$Height_tplus1*coeff_bigger,  pch=1, col=my.palette, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger, c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )

–g1_10: mu_x = b0 + b1(diameter) + b2(height) + epsilon(biotype{E,H,W})

         mu_y = b3 + b4(diameter) + b5(height) + epsilon(biotype{E,H,W})

g1_10 <-manova(cbind(seedlings$Diameter_tplus1, seedlings$Height_tplus1) ~ seedlings$Diameter_t+seedlings$Height_t + seedlings$Genetic_type)
summary(g1_10)

##                         Df   Pillai approx F num Df den Df    Pr(>F)    
## seedlings$Diameter_t     1 0.285801   69.029      2    345 < 2.2e-16 ***
## seedlings$Height_t       1 0.161260   33.166      2    345 6.695e-14 ***
## seedlings$Genetic_type   2 0.050986    4.526      4    692  0.001293 ** 
## Residuals              346                                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sd_x<-summary(lm(seedlings$Diameter_tplus1 ~ seedlings$Diameter_t + seedlings$Height_t + seedlings$Genetic_type))$sigma

##EASTERN
#Grab parameters associated with future diameter
b0<-g1_10$coeff[1,1]
b1<-g1_10$coeff[2,1]
b2<-g1_10$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_10$coeff[1,2]
b4<-g1_10$coeff[2,2]
b5<-g1_10$coeff[3,2]
sd_y <-summary(lm(seedlings$Height_tplus1 ~ seedlings$Diameter_t + seedlings$Height_t + seedlings$Genetic_type))$sigma
 
 

z1_E<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_E<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

par(mfrow=c(1,2))
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_E, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Eastern", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings$Diameter_t[seedlings$Genetic_type == "eastern"], seedlings$Height_t[seedlings$Genetic_type == "eastern"], seedlings$Diameter_tplus1[seedlings$Genetic_type == "eastern"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_E, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Eastern", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings$Diameter_t[seedlings$Genetic_type == "eastern"], seedlings$Height_t[seedlings$Genetic_type == "eastern"], seedlings$Height_tplus1[seedlings$Genetic_type =="eastern"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

##Hybrid
#Grab parameters associated with future diameter
b0<-g1_10$coeff[1,1] + g1_10$coeff[4,1]
b1<-g1_10$coeff[2,1]
b2<-g1_10$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_10$coeff[1,2] + g1_10$coeff[4,2]
b4<-g1_10$coeff[2,2]
b5<-g1_10$coeff[3,2]
sd_y <-summary(lm(seedlings$Height_tplus1 ~ seedlings$Diameter_t + seedlings$Height_t + seedlings$Genetic_type))$sigma
 
 

z1_H<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_H<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

par(mfrow=c(1,2))
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_H, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Hybrid", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings$Diameter_t[seedlings$Genetic_type == "hybrid"], seedlings$Height_t[seedlings$Genetic_type == "hybrid"], seedlings$Diameter_tplus1[seedlings$Genetic_type == "hybrid"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_H, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Hybrid", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings$Diameter_t[seedlings$Genetic_type == "hybrid"], seedlings$Height_t[seedlings$Genetic_type == "hybrid"], seedlings$Height_tplus1[seedlings$Genetic_type =="hybrid"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

###Western
#Grab parameters associated with future diameter
b0<-g1_10$coeff[1,1] + g1_10$coeff[5,1]
b1<-g1_10$coeff[2,1]
b2<-g1_10$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_10$coeff[1,2] + g1_10$coeff[5,2]
b4<-g1_10$coeff[2,2]
b5<-g1_10$coeff[3,2]
sd_y <-summary(lm(seedlings$Height_tplus1 ~ seedlings$Diameter_t + seedlings$Height_t + seedlings$Genetic_type))$sigma
 
 

z1_W<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_W<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

par(mfrow=c(1,2))
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_W, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Western", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings$Diameter_t[seedlings$Genetic_type == "western"], seedlings$Height_t[seedlings$Genetic_type == "western"], seedlings$Diameter_tplus1[seedlings$Genetic_type == "western"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_W, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Western", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings$Diameter_t[seedlings$Genetic_type == "western"], seedlings$Height_t[seedlings$Genetic_type == "western"], seedlings$Height_tplus1[seedlings$Genetic_type =="western"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

#EASTERN IMAGE PLOT
par(mfrow=c(1,2))
my.palette <- viridis(9)
par(mar=c(5.1,4.1,4.1,6))
coeff_bigger=.1
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_E, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "Eastern")
points(seedlings$Diameter_t[seedlings$Genetic_type == "eastern"], seedlings$Height_t[seedlings$Genetic_type == "eastern"], cex=seedlings$Diameter_tplus1[seedlings$Genetic_type == "eastern"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")
#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_E, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="Eastern",legend.lab = "Height at t+1 (cm)")
points(seedlings$Diameter_t[seedlings$Genetic_type == "eastern"], seedlings$Height_t[seedlings$Genetic_type == "eastern"], cex=seedlings$Height_tplus1[seedlings$Genetic_type == "eastern"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )

#Hybrid IMAGE PLOT
par(mfrow=c(1,2))
my.palette <- viridis(9)
par(mar=c(5.1,4.1,4.1,6))
coeff_bigger=.1
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_H, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "Hybrid")
points(seedlings$Diameter_t[seedlings$Genetic_type == "hybrid"], seedlings$Height_t[seedlings$Genetic_type == "hybrid"], cex=seedlings$Diameter_tplus1[seedlings$Genetic_type == "hybrid"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")
#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_H, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="Hybrid",legend.lab = "Height at t+1 (cm)")
points(seedlings$Diameter_t[seedlings$Genetic_type == "hybrid"], seedlings$Height_t[seedlings$Genetic_type == "hybrid"], cex=seedlings$Height_tplus1[seedlings$Genetic_type == "hybrid"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )

#Western IMAGE PLOT
par(mfrow=c(1,2))
my.palette <- viridis(9)
par(mar=c(5.1,4.1,4.1,6))
coeff_bigger=.1
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_W, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "Western")
points(seedlings$Diameter_t[seedlings$Genetic_type == "western"], seedlings$Height_t[seedlings$Genetic_type == "western"], cex=seedlings$Diameter_tplus1[seedlings$Genetic_type == "western"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")
#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_W, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="Western",legend.lab = "Height at t+1 (cm)")
points(seedlings$Diameter_t[seedlings$Genetic_type == "western"], seedlings$Height_t[seedlings$Genetic_type == "western"], cex=seedlings$Height_tplus1[seedlings$Genetic_type == "western"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )

–g1_11: mu_x = b0 + b1(diameter) + b2(height) + epsilon(site{BC, CC, C, FP, PG,WT})

         mu_y = b0 + b1(diameter) + b2(height) + epsilon(site{BC, CC, C, FP, PG, WT})

g1_11 <-manova(cbind(seedlings_g$Diameter_tplus1, seedlings_g$Height_tplus1) ~ seedlings_g$Diameter_t+seedlings_g$Height_t + seedlings_g$Site)
summary(g1_11)

##                         Df  Pillai approx F num Df den Df    Pr(>F)    
## seedlings_g$Diameter_t   1 0.33573   59.640      2    236 < 2.2e-16 ***
## seedlings_g$Height_t     1 0.09784   12.797      2    236 5.292e-06 ***
## seedlings_g$Site         5 0.12281    3.101     10    474 0.0007681 ***
## Residuals              237                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sd_x<-summary(lm(seedlings_g$Diameter_tplus1 ~ seedlings_g$Diameter_t + seedlings_g$Height_t + seedlings_g$Site))$sigma

##Big Cypress
#Grab parameters associated with future diameter
b0<-g1_11$coeff[1,1]
b1<-g1_11$coeff[2,1]
b2<-g1_11$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_11$coeff[1,2]
b4<-g1_11$coeff[2,2]
b5<-g1_11$coeff[3,2]
sd_y <-summary(lm(seedlings$Height_tplus1 ~ seedlings$Diameter_t + seedlings$Height_t + seedlings$Site))$sigma
 
 

z1_BC<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_BC<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

##Cape Canaveral
#Grab parameters associated with future diameter
b0<-g1_11$coeff[1,1] + g1_11$coeff[4,1]
b1<-g1_11$coeff[2,1]
b2<-g1_11$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_11$coeff[1,2] + g1_11$coeff[4,2]
b4<-g1_11$coeff[2,2]
b5<-g1_11$coeff[3,2]

z1_CC<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_CC<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

##Chekika
#Grab parameters associated with future diameter
b0<-g1_11$coeff[1,1] + g1_11$coeff[5,1]
b1<-g1_11$coeff[2,1]
b2<-g1_11$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_11$coeff[1,2] + g1_11$coeff[5,2]
b4<-g1_11$coeff[2,2]
b5<-g1_11$coeff[3,2]

z1_C<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_C<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

##Fort Pierce
#Grab parameters associated with future diameter
b0<-g1_11$coeff[1,1] + g1_11$coeff[6,1]
b1<-g1_11$coeff[2,1]
b2<-g1_11$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_11$coeff[1,2] + g1_11$coeff[6,2]
b4<-g1_11$coeff[2,2]
b5<-g1_11$coeff[3,2]

z1_FP<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_FP<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))


##Punta Gorda
#Grab parameters associated with future diameter
b0<-g1_11$coeff[1,1] + g1_11$coeff[7,1]
b1<-g1_11$coeff[2,1]
b2<-g1_11$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_11$coeff[1,2] + g1_11$coeff[7,2]
b4<-g1_11$coeff[2,2]
b5<-g1_11$coeff[3,2]

z1_PG<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_PG<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

##Wild Turkey
#Grab parameters associated with future diameter
b0<-g1_11$coeff[1,1] + g1_11$coeff[8,1]
b1<-g1_11$coeff[2,1]
b2<-g1_11$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_11$coeff[1,2] + g1_11$coeff[8,2]
b4<-g1_11$coeff[2,2]
b5<-g1_11$coeff[3,2]

z1_WT<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_WT<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))


##### 3D plots
par(mfrow=c(1, 2))

###Big Cypress
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_BC, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Big Cypress", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Big Cypress"], seedlings_g$Height_t[seedlings_g$Site == "Big Cypress"], seedlings_g$Diameter_tplus1[seedlings_g$Site == "Big Cypress"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_BC, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Big Cypress", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Big Cypress"], seedlings_g$Height_t[seedlings_g$Site == "Big Cypress"], seedlings_g$Height_tplus1[seedlings_g$Site =="Big Cypress"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

###Cape Canaveral
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_CC, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Cape Canaveral", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Cape Canaveral"], seedlings_g$Height_t[seedlings$Site == "Cape Canaveral"], seedlings_g$Diameter_tplus1[seedlings$Site == "Cape Canaveral"], pmat=pmat)

## Warning in cbind(x, y, z, 1, deparse.level = 0L): number of rows of result
## is not a multiple of vector length (arg 1)

# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_CC, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Cape Canaveral", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings$Diameter_t[seedlings$Site == "Cape Canaveral"], seedlings$Height_t[seedlings$Site == "Cape Canaveral"], seedlings$Height_tplus1[seedlings$Site =="Cape Canaveral"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

###Chekika
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_C, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Chekika", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Chekika"], seedlings_g$Height_t[seedlings_g$Site == "Chekika"], seedlings_g$Diameter_tplus1[seedlings_g$Site == "Chekika"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_C, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Chekika", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Chekika"], seedlings_g$Height_t[seedlings_g$Site == "Chekika"], seedlings_g$Height_tplus1[seedlings_g$Site =="Chekika"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

###Fort Pierce
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_FP, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Fort Pierce", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Fort Pierce"], seedlings_g$Height_t[seedlings_g$Site == "Fort Pierce"], seedlings_g$Diameter_tplus1[seedlings_g$Site == "Fort Pierce"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_FP, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Fort Pierce", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Fort Pierce"], seedlings_g$Height_t[seedlings_g$Site == "Fort Pierce"], seedlings_g$Height_tplus1[seedlings_g$Site =="Fort Pierce"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

### Punta Gorda
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_PG, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Punta Gorda", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Punta Gorda"], seedlings_g$Height_t[seedlings_g$Site == "Punta Gorda"], seedlings_g$Diameter_tplus1[seedlings_g$Site == "Punta Gorda"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_PG, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Punta Gorda", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Punta Gorda"], seedlings_g$Height_t[seedlings_g$Site == "Punta Gorda"], seedlings_g$Height_tplus1[seedlings_g$Site =="Punta Gorda"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

### Wild Turkey
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_WT, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Wild Turkey", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Wild Turkey"], seedlings_g$Height_t[seedlings_g$Site == "Wild Turkey"], seedlings_g$Diameter_tplus1[seedlings_g$Site == "Wild Turkey"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_WT, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="Wild Turkey", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings_g$Diameter_t[seedlings_g$Site == "Wild Turkey"], seedlings_g$Height_t[seedlings_g$Site == "Wild Turkey"], seedlings_g$Height_tplus1[seedlings_g$Site =="Wild Turkey"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

##### Image Plots
par(mfrow=c(1,2))
my.palette <- viridis(9)
par(mar=c(5.1,4.1,4.1,6))
coeff_bigger=.1

#Big Cypress IMAGE PLOT
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_BC, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "Big Cypress")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Big Cypress"], seedlings_g$Height_t[seedlings_g$Site == "Big Cypress"], cex=seedlings_g$Diameter_tplus1[seedlings_g$Site == "Big Cypress"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_BC, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="Big Cypress",legend.lab = "Height at t+1 (cm)")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Big Cypress"], seedlings_g$Height_t[seedlings_g$Site == "Big Cypress"], cex=seedlings_g$Height_tplus1[seedlings_g$Site == "Big Cypress"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#Cape Canaveral IMAGE PLOT
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_CC, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "Cape Canaveral")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Cape Canaveral"], seedlings_g$Height_t[seedlings_g$Site == "Cape Canaveral"], cex=seedlings_g$Diameter_tplus1[seedlings_g$Site == "Cape Canaveral"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_CC, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="Cape Canaveral",legend.lab = "Height at t+1 (cm)")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Cape Canaveral"], seedlings_g$Height_t[seedlings_g$Site == "Cape Canaveral"], cex=seedlings_g$Height_tplus1[seedlings_g$Site == "Cape Canaveral"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#Chekika IMAGE PLOT
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_C, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "Chekika")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Chekika"], seedlings_g$Height_t[seedlings_g$Site == "Chekika"], cex=seedlings_g$Diameter_tplus1[seedlings_g$Site == "Chekika"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_C, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="Chekika",legend.lab = "Height at t+1 (cm)")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Chekika"], seedlings_g$Height_t[seedlings_g$Site == "Chekika"], cex=seedlings_g$Height_tplus1[seedlings_g$Site == "Chekika"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#Fort Pierce IMAGE PLOT
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_FP, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "Fort Pierce")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Fort Pierce"], seedlings_g$Height_t[seedlings_g$Site == "Fort Pierce"], cex=seedlings_g$Diameter_tplus1[seedlings_g$Site == "Fort Pierce"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_FP, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="Fort Pierce",legend.lab = "Height at t+1 (cm)")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Fort Pierce"], seedlings_g$Height_t[seedlings_g$Site == "Fort Pierce"], cex=seedlings_g$Height_tplus1[seedlings_g$Site == "Fort Pierce"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#Punta Gorda IMAGE PLOT
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_PG, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "Punta Gorda")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Punta Gorda"], seedlings_g$Height_t[seedlings_g$Site == "Punta Gorda"], cex=seedlings_g$Diameter_tplus1[seedlings_g$Site == "Punta Gorda"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_PG, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="Punta Gorda",legend.lab = "Height at t+1 (cm)")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Punta Gorda"], seedlings_g$Height_t[seedlings_g$Site == "Punta Gorda"], cex=seedlings_g$Height_tplus1[seedlings_g$Site == "Punta Gorda"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#Wild Turkey IMAGE PLOT
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_WT, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "Wild Turkey")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Wild Turkey"], seedlings_g$Height_t[seedlings_g$Site == "Wild Turkey"], cex=seedlings_g$Diameter_tplus1[seedlings_g$Site == "Wild Turkey"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_WT, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="Wild Turkey",legend.lab = "Height at t+1 (cm)")
points(seedlings_g$Diameter_t[seedlings_g$Site == "Wild Turkey"], seedlings_g$Height_t[seedlings_g$Site == "Wild Turkey"], cex=seedlings_g$Height_tplus1[seedlings_g$Site == "Wild Turkey"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

–g1_12: mu_x = b0 + b1(diameter) + b2(height) + epsilon(year{’10, ’11, ’12, ’13, ’14})

         mu_y = b3 + b4(diameter) + b5(height) + epsilon(year{'10, '11, '12, '13, '14})

g1_12 <-manova(cbind(seedlings2$Diameter_tplus1, seedlings2$Height_tplus1) ~ seedlings2$Diameter_t+seedlings2$Height_t + seedlings2$Year)
summary(g1_12)

##                        Df  Pillai approx F num Df den Df    Pr(>F)    
## seedlings2$Diameter_t   1 0.21563   33.676      2    245 1.200e-13 ***
## seedlings2$Height_t     1 0.16754   24.654      2    245 1.756e-10 ***
## seedlings2$Year         3 0.15517    6.897      6    492 4.780e-07 ***
## Residuals             246                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sd_x<-summary(lm(seedlings2$Diameter_tplus1 ~ seedlings2$Diameter_t + seedlings2$Height_t + seedlings2$Year))$sigma
sd_y <-summary(lm(seedlings2$Height_tplus1 ~ seedlings2$Diameter_t + seedlings2$Height_t + seedlings2$Year))$sigma

##2010
#Grab parameters associated with future diameter
b0<-g1_12$coeff[1,1]
b1<-g1_12$coeff[2,1]
b2<-g1_12$coeff[3,1]

#Grab parameters associated with future height
b3<-g1_12$coeff[1,2]
b4<-g1_12$coeff[2,2]
b5<-g1_12$coeff[3,2]

z1_2010<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_2010<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

##2011
#Grab parameters associated with future diameter
b0<-g1_12$coeff[1,1] + g1_12$coeff[4,1]
b1<-g1_12$coeff[2,1]
b2<-g1_12$coeff[3,1]

#Grab parameters associated with future height
b3<-g1_12$coeff[1,2] + g1_12$coeff[4,2]
b4<-g1_12$coeff[2,2]
b5<-g1_12$coeff[3,2]

z1_2011<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_2011<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

##2012
#Grab parameters associated with future diameter
b0<-g1_12$coeff[1,1] + g1_12$coeff[5,1]
b1<-g1_12$coeff[2,1]
b2<-g1_12$coeff[3,1]

#Grab parameters associated with future height
b3<-g1_12$coeff[1,2] + g1_12$coeff[5,2]
b4<-g1_12$coeff[2,2]
b5<-g1_12$coeff[3,2]

z1_2012<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_2012<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

##2013
#Grab parameters associated with future diameter
b0<-g1_12$coeff[1,1] + g1_12$coeff[6,1]
b1<-g1_12$coeff[2,1]
b2<-g1_12$coeff[3,1]

#Grab parameters associated with future height
b3<-g1_12$coeff[1,2] + g1_12$coeff[6,2]
b4<-g1_12$coeff[2,2]
b5<-g1_12$coeff[3,2]

z1_2013<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b )))
z2_2013<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

par(mfrow=c(1,2))
#3D Perspective Plots 
### 2010
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_2010, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="2010-2011", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings2$Diameter_t[seedlings2$Year == "2010"], seedlings2$Height_t[seedlings2$Year == "2010"], seedlings2$Diameter_tplus1[seedlings2$Year == "2010"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_2010, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="2010-2011", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings2$Diameter_t[seedlings2$Year == "2010"], seedlings2$Height_t[seedlings2$Year == "2010"], seedlings2$Height_tplus1[seedlings2$Year =="2010"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

### 2011
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_2011, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="2011-2012", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings2$Diameter_t[seedlings2$Year == "2011"], seedlings2$Height_t[seedlings2$Year == "2011"], seedlings2$Diameter_tplus1[seedlings2$Year == "2011"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_2011, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="2011-2012", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings2$Diameter_t[seedlings2$Year == "2011"], seedlings2$Height_t[seedlings2$Year == "2011"], seedlings2$Height_tplus1[seedlings2$Year =="2011"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

### 2012
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_2012, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="2012-2013", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings2$Diameter_t[seedlings2$Year == "2012"], seedlings2$Height_t[seedlings2$Year == "2012"], seedlings2$Diameter_tplus1[seedlings2$Year == "2012"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_2012, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="2012-2013", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings2$Diameter_t[seedlings2$Year == "2012"], seedlings2$Height_t[seedlings2$Year == "2012"], seedlings2$Height_tplus1[seedlings2$Year =="2012"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

### 2013
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1_2013, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="2013-2014", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings2$Diameter_t[seedlings2$Year == "2013"], seedlings2$Height_t[seedlings2$Year == "2013"], seedlings2$Diameter_tplus1[seedlings2$Year == "2013"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)


#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2_2013, ticktype="detailed",
      theta=-30, zlim=c(0,44), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="2013-2014", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings2$Diameter_t[seedlings2$Year == "2013"], seedlings2$Height_t[seedlings2$Year == "2013"], seedlings2$Height_tplus1[seedlings2$Year =="2013"], pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

##### Image Plots 
#2010
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_2010, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "2010")
points(seedlings2$Diameter_t[seedlings2$Year == "2010"], seedlings2$Height_t[seedlings2$Year == "2010"], cex=seedlings2$Diameter_tplus1[seedlings2$Year == "2010"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_2010, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="2010",legend.lab = "Height at t+1 (cm)")
points(seedlings2$Diameter_t[seedlings2$Year == "2010"], seedlings2$Height_t[seedlings2$Year == "2010"], cex=seedlings2$Height_tplus1[seedlings2$Year == "2010"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#2011
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_2011, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "2011")
points(seedlings2$Diameter_t[seedlings2$Year == "2011"], seedlings2$Height_t[seedlings2$Year == "2011"], cex=seedlings2$Diameter_tplus1[seedlings2$Year == "2011"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_2011, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="2011",legend.lab = "Height at t+1 (cm)")
points(seedlings2$Diameter_t[seedlings2$Year == "2011"], seedlings2$Height_t[seedlings2$Year == "2011"], cex=seedlings2$Height_tplus1[seedlings2$Year == "2011"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#2012
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_2012, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "2012")
points(seedlings2$Diameter_t[seedlings2$Year == "2012"], seedlings2$Height_t[seedlings2$Year == "2012"], cex=seedlings2$Diameter_tplus1[seedlings2$Year == "2012"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_2012, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="2012",legend.lab = "Height at t+1 (cm)")
points(seedlings2$Diameter_t[seedlings2$Year == "2012"], seedlings2$Height_t[seedlings2$Year == "2012"], cex=seedlings2$Height_tplus1[seedlings2$Year == "2012"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#2013
#Future Diameter
image.plot(x1seq_seedlings, x2seq_seedlings, z1_2013, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)", main = "2013")
points(seedlings2$Diameter_t[seedlings2$Year == "2013"], seedlings2$Height_t[seedlings2$Year == "2013"], cex=seedlings2$Diameter_tplus1[seedlings2$Year == "2013"]*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#Future Height
image.plot(x1seq_seedlings, x2seq_seedlings, z2_2013, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, main ="2013",legend.lab = "Height at t+1 (cm)")
points(seedlings2$Diameter_t[seedlings2$Year == "2013"], seedlings2$Height_t[seedlings2$Year == "2013"], cex=seedlings2$Height_tplus1[seedlings2$Year == "2013"]*coeff_bigger,  pch=1, bg = "white")
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger,c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

- AIC

#AIC is not very useful here, since the different classes of models are modeling different things. And AIC cannot be computed for the models where both diameter and height are being modeled?

#For all models including diameter only (no year)
aictab(list(g1_1, g1_2, g1_3), modnames = c("g1_1", "g1_2", "g1_3"))

#For all models including height only
aictab(list(g1_5, g1_6, g1_7), modnames = c("g1_5", "g1_6", "g1_7"))

2. Corrected to only include individuals within D1

–g1_1c Diameter only

seedlings3 <- subset(seedlings, seedlings$Diameter_tplus1 < 1.6 & seedlings$Height_tplus1 < 16)
seedlings4 <- subset(seedlings, seedlings$Diameter_tplus1 >1.6 & seedlings$Height_tplus1>16)

#Seedlings that remain in the seedling domain 
g1_1c <- lm(seedlings3$Diameter_tplus1 ~ seedlings3$Diameter_t)
summary(g1_1c)

## 
## Call:
## lm(formula = seedlings3$Diameter_tplus1 ~ seedlings3$Diameter_t)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83841 -0.12663 -0.01484  0.10569  0.54105 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.56769    0.04895   11.60   <2e-16 ***
## seedlings3$Diameter_t  0.55894    0.05317   10.51   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1927 on 243 degrees of freedom
## Multiple R-squared:  0.3126, Adjusted R-squared:  0.3098 
## F-statistic: 110.5 on 1 and 243 DF,  p-value: < 2.2e-16

#Variance in Growth
sd_x <- summary(g1_1c)$sigma


preds2 <- g1_1c$coefficients[1] + g1_1c$coefficients[2]*x1seq_seedlings + 2*sd_x
preds3 <- g1_1c$coefficients[1] + g1_1c$coefficients[2]*x1seq_seedlings - 2*sd_x
preds4 <- g1_1c$coefficients[1] + g1_1c$coefficients[2]*x1seq_seedlings + 1*sd_x
preds5 <- g1_1c$coefficients[1] + g1_1c$coefficients[2]*x1seq_seedlings - 1*sd_x



plot(x1seq_seedlings, g1_1c$coeff[1] + g1_1c$coeff[2]*x1seq_seedlings, type='l' , xlab = "Diameter at t (mm)", ylab="Diameter at t+1 (mm)", main = "Seedlings that remain seedlings")
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x1seq_seedlings), x1seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings3$Diameter_t, seedlings3$Diameter_tplus1)
abline(g1_1c)

–g1_5c Height only

#Seedlings that remain in the seedling domain 
g1_5c <- lm(seedlings3$Height_tplus1 ~ seedlings3$Height_t)
summary(g1_5c)

## 
## Call:
## lm(formula = seedlings3$Height_tplus1 ~ seedlings3$Height_t)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.654 -1.887  0.075  1.729  5.767 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          6.81236    0.60623  11.237  < 2e-16 ***
## seedlings3$Height_t  0.34585    0.06688   5.171 4.87e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.713 on 243 degrees of freedom
## Multiple R-squared:  0.09913,    Adjusted R-squared:  0.09542 
## F-statistic: 26.74 on 1 and 243 DF,  p-value: 4.868e-07

#Variance in Growth
sd_x <- summary(g1_5c)$sigma


preds2 <- g1_5c$coefficients[1] + g1_5c$coefficients[2]*x2seq_seedlings + 2*sd_x
preds3 <- g1_5c$coefficients[1] + g1_5c$coefficients[2]*x2seq_seedlings - 2*sd_x
preds4 <- g1_5c$coefficients[1] + g1_5c$coefficients[2]*x2seq_seedlings + 1*sd_x
preds5 <- g1_5c$coefficients[1] + g1_5c$coefficients[2]*x2seq_seedlings - 1*sd_x



plot(x2seq_seedlings, g1_5c$coeff[1] + g1_5c$coeff[2]*x2seq_seedlings, type='l' , xlab = "Height at t (cm)", ylab="Height at t+1 (cm)", main = "Seedlings that remain seedlings")
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds3), preds2), col = 'grey90', border = NA)
polygon(c(rev(x2seq_seedlings), x2seq_seedlings), c(rev(preds5), preds4), col = 'grey80', border = NA)
points(seedlings3$Height_t, seedlings3$Height_tplus1)
abline(g1_5c)

#–g1_9c Diameter + Height #This is still not doing a great job (notice scale). The model is vastly underestimating growth.

  g1_9c <-manova(cbind(seedlings3$Diameter_tplus1, seedlings3$Height_tplus1) ~ seedlings3$Diameter_t+seedlings3$Height_t)
summary(g1_9c)

##                        Df  Pillai approx F num Df den Df    Pr(>F)    
## seedlings3$Diameter_t   1 0.31903   56.452      2    241 < 2.2e-16 ***
## seedlings3$Height_t     1 0.09554   12.729      2    241 5.558e-06 ***
## Residuals             242                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sd_x<-summary(lm(seedlings3$Diameter_tplus1 ~ seedlings3$Diameter_t + seedlings3$Height_t))$sigma

#Grab parameters associated with future diameter
b0<-g1_9c$coeff[1,1]
b1<-g1_9c$coeff[2,1]
b2<-g1_9c$coeff[3,1]



#Grab parameters associated with future height
b3<-g1_9c$coeff[1,2]
b4<-g1_9c$coeff[2,2]
b5<-g1_9c$coeff[3,2]
sd_y <-summary(lm(seedlings3$Height_tplus1 ~ seedlings3$Diameter_t + seedlings3$Height_t))$sigma
 
 

z1<-outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b0 + b1*a + b2*b)))
z2<- outer(x1seq_seedlings, x2seq_seedlings, function(a,b) (invlogit(b3+ b4*a + b5*b)))

par(mfrow=c(1,2))
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z1, ticktype="detailed",
      theta=-30, zlim=c(0,4.1), zlab="\n Diameter at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings3$Diameter_t, seedlings3$Height_t, seedlings3$Diameter_tplus1, pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)
#Now for future height:
pmat<-persp(x1seq_seedlings, x2seq_seedlings, z2, ticktype="detailed",
      theta=-30, zlim=c(0,16), zlab="\n Height at t+1 (mm)", 
      shade=0.1, nticks=4, xlab="\n Diameter at t(mm)", 
      ylab="\n Height at t(cm)",  lwd=1, xaxs="i", 
      main="", cex.main=1)

# from 3D to 2D coordinates
mypoints <- trans3d(seedlings3$Diameter_t, seedlings3$Height_t, seedlings3$Height_tplus1, pmat=pmat)
# plot in 2D space with pointsize related to distance
points(mypoints, col=4)

par(mfrow=c(1,1))
par(mar=c(5.1,4.1,4.1,6))
coeff_bigger=1
image.plot(x1seq_seedlings, x2seq_seedlings, z1, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Diameter at t+1 (mm)")
points(seedlings3$Diameter_t, seedlings3$Height_t, cex=seedlings3$Diameter_tplus1*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(0.5,1, 1.5, 2,2.5, 3, 3.5,4), c("0.5", "1", "1.5", "2", "2.5", "3", "3.5", "4"), xpd=T, horiz = T, title="Diameter at t +1 (mm)")

#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )


my.palette <- viridis(9)
par(mar=c(5.1,4.1,4.1,6))
coeff_bigger=0.1
image.plot(x1seq_seedlings, x2seq_seedlings, z2, xlab = "Diameter (mm) ", ylab="Height (cm)", col=my.palette, legend.lab = "Height at t+1 (cm)")
points(seedlings3$Diameter_t, seedlings3$Height_t, cex=seedlings3$Height_tplus1*coeff_bigger,  pch=1)
legend(0, 18.5, pch=1,text.width=0.1, pt.cex = c(3,8, 13,18, 23, 28, 33, 38, 44)*coeff_bigger, c("3", "8", "13", "23", "28", "33", "38", "44"), xpd=T, horiz = T, title="Height at t +1 (cm)")

#text (seedlings$Diameter_t , seedlings$Height_t , seedlings$Diameter_tplus1 )

Exploring vital rate models: D1 Growth

Kelley Erickson

May 16, 2019